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There is increasing interest in the analysis of biological tissue, its organization and its 
dynamics with the help of mathematical models. In the ideal case emergent properties on 
the tissue scale can be derived from the cellular scale. However, this has been achieved 
in rare examples only, in particular, when involving high-speed migration of cells. One 
major difficulty is the lack of a suitable multiscale simulation platform, which embeds 
reaction-diffusion of soluble substances, fast cell migration and mechanics, and, being of 
great importance in several tissue types, cell flow homeostasis. In this paper a step into 
this direction is presented by developing an agent-based mathematical model specifically 
designed to incorporate these features with special emphasis on high speed cell migration. 
Cells are represented as elastic spheres migrating on a substrate in lattice-free space. Their 
movement is regulated and guided by chemoattractants that can be derived from the 
substrate. The diffusion of chemoattractants is considered to be slower than cell migration 
and, thus, to be far from equilibrium. Tissue homeostasis is not achieved by the balance 
of growth and death but by a flow equilibrium of cells migrating in and out of the tissue 
under consideration. In this sense the number and the distribution of the cells in the tissue 
is a result of the model and not part of the assumptions. For purpose of demonstration 
of the model properties and functioning, the model is applied to a prominent example of 
tissue in a cellular flow equilibrium, the secondary lymphoid tissue. The experimental data 
on cell speed distributions in these tissues can be reproduced using reasonable mechanical 
parameters for the simulated cell migration in dense tissue. 
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1. Introduction 

There exists a number of agent-based off-lattice models to simulate tissue [1-8]. Many of them describe 
tissue organization and pattern formation [9] by the dynamics of growth and death of cells. Other 
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models are suitable for describing migration of a constant number of cells. The traditional research field 
of pattern formation involving migrating cells [3,4,10] is mostly based on a balance of proliferation and 
cell death. Cells are migrating with relatively slow speed with respect to other processes like diffusion 
of soluble substances. Approximations based on slow cell migration are justified even for such highly 
dynamic systems as the epidermis [8, 11]. 

The model introduced here is designed for tissue in a flow equilibrium of high-speed migrating 
cells as found frequently in immunological tissue. Cells can enter and leave the tissue according 
to some dynamics imposed by the considered biological system. Feedback mechanisms between the 
cells and the substrate are incorporated because homeostasis in such systems is often established by 
feedback interaction [12-14]. This is modeled as transient excitation of a substrate or a small sessile 
cell population. The excitation is thought to couple back to the cell migration thus allowing for a 
regulatory system in a stable flow equilibrium. 

The model is based on a regular triangulation [15, 16] that provides the cell neighborship topology 
which is changing rapidly due to fast cell migration. The dual Voronoi tessellation [16,17] provides in- 
formation about the shape of the cells including cell contact surfaces [8,11,18]. The physical properties 
are based on previously introduced short-ranged elastic interactions and actively generated intercel- 
lular forces [3-5,8,18]. Each cell is modeled individually and can incorporate a set of molecular 
properties or internal states that is appropriate for the system under consideration. Such properties 
may have impact on cell migration or cell-cell and cell-substrate interaction. The deterministic internal 
cell dynamics exhibit memory in the sense that, for example, internal variables representing the cell 
state like delay times are included. Cell migration under the influence of chemotaxis is also treated 
deterministic. However, cell orientation is stochastic for unguided random cell migration. 

In order to illustrate the model features, the model is used to investigate certain features of cell 
behavior that are relevant for the formation of primary lymphoid follicles (PLF). These can be found 
in secondary lymphoid tissue (SLT) like the spleen and lymph nodes [12,14,19-22]. Two lymphocyte 
cell populations of the PLF show fast migration: B cells and T cells. Two sessile populations, namely 
follicular dendritic cells (FDC) and fibroblastic reticular cells (FRC) serve as a complex 3D substrate 
for lymphocytes. The morphology of SLT is characterized by PLF and T zones. The PLF contains the 
B cells and the FDC. It is adjacent to the T zone harboring T cells and FRC. The functionality of the 
model features is demonstrated by analyzing the migration of the lymphocytes and their interaction 
with chemokines. Also the generation of a tissue in flow equilibrium is studied using a toy model 
inspired by PLF formation. On the basis of biologically proven interactions the model can explain a 
clear separation of T and B cells forming a T zone and an FDC containing follicle, as found in real 
tissue. However, more realistic models of PLF remain to be developed in future work. 

The presented simulation technique has a direct connection to experimental constraints and exhibits 
the potential of building physically concise models of cell migration in tissue. It's novelty with respect 
to previous work relies on the focus on tissues composed of rapidly migrating cells which exhibit 
cell-flow-equilibrium. This extends the range of applications of agent-based off-lattice techniques to 
lymphoid tissue. 
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2. Method 



The tissue formation problem of migrating cells is simulated using an agent-based model on top of 
a regular triangulation [15,16,23,24]. The regular triangulation is used to provide the neighborhood 
topology for the cells that allows for a continuous representation of cell positions and sizes in contrast 
to grid-based methods. The regular triangulation also provides information about the cell contacts 
and contact areas, and, within the limits of an approximation, about cell volume and shape [8,25]. 
The regular triangulation will be briefly summarized in the next subsection. 

The simulation of cells is realized in a 3-level multiscale model. The first level is the internal state 
of the cells representing the dynamics of the phenotype of the cell (see Sec. I2.2p . The second level 
models the contact interaction between cells including mechanical interactions with the environment 
and exchange of signals by membrane bound molecules (see Sec. 12. 3D . The third level incorporates 
long range interactions via diffusive substances for example chemotaxis, i.e. directed motion of cells 
upwards a concentration gradient which is induced by molecular chemoattractants (see Sec. l2.4h . Thus, 
chemoattractants are derived from cellular sources and feed back to cell migration. 



2.1. Regular triangulation 

Each cell is represented as a sphere at position x with radius R. A vertex is defined by the pair 
X = (x,R). The regular triangulation is defined using the empty orthosphere criterion [15-17]. 
In three dimensions four vertices A,B,C,D forming a tetrahedron uniquely define an orthosphere 
(Fig. [I]). The orthosphere is empty if for any other vertex V 
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holds provided that the four vertices are oriented positively, i.e. if 
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The symbols a x 



denote the coordinates of vertex A at the position a 



,a y ,a z ) and R a the 



corresponding radius of the elastic sphere associated with the vertex. The notation for the vertices 
B,C,D,V is analogously defined. A set of non-overlapping tetrahedras covering a set of vertices 
forms a regular triangulation if all orthospheres attributed to these tetrahedras do not contain any 
further vertex. The regular triangulation, and with it the neighborhood topology, changes when cells 
are moving, and when they are added or removed from the system due to cell flux, cell death or 
cell proliferation. The corresponding algorithms have been developed and published previously for 
serial [23] and parallel computer architectures [24]. 
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Figure 1: Orthosphere in two dimensions (left panel). The triangle ABC defines an orthocircle M 
in two dimensions. The weights of the vertices are represented as disks with corresponding 
radii. Geometrically, the orthogonality of the orthocircle is indicated by the orthogonality 
of the tangents at the intersection point of the orthocircle and the circles representing the 
vertex weight. A vertex V is inside the orthocircle M when the tangential intersection lies 
within M, as shown in the example. A regular triangulation is marked as dashed straight 
lines connecting the vertices in the right panel. The dark grey polygon is the Voronoi cell of 
the central vertex with an associated weighted sphere (grey disk). The dashed circles depict 
the weighted spheres of the neighbor vertices. 

2.2. Internal cell dynamics 

The phenotype of a cell is described by a set of internal cell variables <j). These include internal times 
to indicate when which type of event may happen. An example is the persistence time T p of cells 
during chemotactic motion [26-30]. When the persistence time T p has past, the cell can reorient to 
the local chemoattractant field [31,32]. 

Upon contact two cells can exchange signals via the contact surface. In a simple approach with 
unpolarized ligand/receptor distributions the signal strength is proportional to the contact area. The 
contact area is computed as the minimum between a sphere overlap and the common Voronoi face 
of these cells (Fig. [2]). This choice relies on the fact that the Voronoi-contact area is the better 
description for high density of cells while the virtual overlap of the spheres is more realistic for low 
density systems [8,18]. In both cases the alternative measure for the shape leads to larger estimates 
which justifies the present choice. In general the minimum of the two contact areas is close to the 
realistic contact of two elastically deforming spherical cells. 

Finally, a set of variables describes the mechanics of the cell: velocity, orientation, cell volume, 
and elasticity. These variable couple directly to the next level of description, the contact interaction 
of cells, and have no direct influence on other internal states but may influence the corresponding 
variables of neighboring cells. 

2.3. Equations of motion 

The contact interaction of cells is predominantly given by mechanical interactions. It is described 
by Newtonian equations of motion in the overdamped approximation [2,4,8]. In this approximation 
acceleration of cells and consequently conservation of moment can be ignored. With dots denoting 
time derivatives we then have 
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Figure 2: Two-dimensional scheme of the contact surface of two cells. The upper panel shows the 
case when the cells are loosely packed. The straight lines indicate the relevant part of the 
Voronoi tessellation. The contact surface between the cells i and j is properly described by 
the overlap a of the two circles (the overlap of two spheres in three dimensions is a disk). 
The lower panel illustrates the case of dense cell packing in which the Voronoi face v is a 
better approximation for the contact surface of the cells i and j. Note, that this definition of 
the contact surface does not only depend on the relative distance of the cells % and j but also 
on cells that are common neighbors of these two cells (indicated by dashed circular lines). 
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The active forces F act on cell i (if any) at position Xj depend on the internal state of the cell and 
the internal state of neighboring cells <j)j, while the passive forces F pass depend on the positions of the 
cell i and its neighbor cells j G Mf in contact with cell i. Active forces can also be exerted directly 
to the surrounding medium. Therefore the set Mi includes all Delaunay-neighbors j independent of 
whether they are in physical contact with the cell i or not. All these forces are counter-balanced by 
the velocity-dependent drag forces F drag . Because of the overdamped approximation this results in an 
ODE system of first order for the cell positions. 

The passive forces are composed of forces stemming from the cell's elasticity and compressibility 

F^ s (x,,x J ) = Ff R (x,,x,)+F— ss ( Xl ,x,) , (4) 
where F^ KR recollects elastic interactions and surface energy. 



2.3.1. Elastic cell-cell interaction 

Elastic forces between cells are treated according to the JKR- model [33]. They depend on the virtual 
cell overlap hij = Ri + Rj — ||xj — Xj|| where Ri and Rj are the cell radii. 
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with the cell's elastic moduli Ei and Ej, Poisson numbers i/j and uj, and the surface energy G{j. The 
force acts in direction of the normal on the contact face 



Ff R (x il x,-)=Jif l (x j ,x i )e 



'■j 



(6) 



The applicability and the limits of the JKR-model for modeling elastic interactions between spherical 
cells have been quantified recently [34]. 



2.3.2. Many-body interactions 

The JKR-model defines a two cell interaction which exhibits an equilibrium distance mediated by 
the balance of surface energy and elastic repulsion. For two cells this leads to a negligible deviation 
of the volume attributed to the vertex and the relaxed volume of the real cell. However, cells will 
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frequently interact with several cells such that specific situations can occur in which cells are strongly 
compressed by surrounding cells without a corresponding large relaxing force being generated. Thus, 
cells might remain in a highly compressed state for too long times because of the neglect of many- 
body interactions. To account for the cell volume and approximately ensure volume conservation, a 
cell pressure concept is included [4, 8] . The pressure of cell i is calculated as deviation of the actual 
cell volume V{ from the target volume V* 

Pi = K \ l ~V*) ( ? ) 
Ki = of , \ s (8) 



3(1 - 2u 



where a linear compression model with compressibility Ki is used. The forces resulting from this 
pressure are exerted between cells by adding the term 

pcompress _ ^ ^ _ ^ 

to the passive cell forces F? ass (xj, xj) in (jlj). is the contact surface of the cells with neighbor cells or 
the medium. In case of contact surface with the medium is calculated as minimum of the Voronoi 
face and the sphere overlap with an equal-sized virtual cell being at equilibrium distance with respect 
to the JKR forces. To get the actual volume the minimum of the sphere volume and the Voronoi cell 
volume is taken which is motivated by similar arguments as for the contact area of two cells [8, 18]. 



2.3.3. Force generation of migrating cells 

Previous models that incorporate the active forces of migrating cells [2-4] use a net propulsion force 
that is exerted into the direction of motion by pulling on neighbor cells. However, cells also exert forces 
perpendicular to the direction of motion which can be even larger than the net force as it has been 
shown for keratocytes [35,36]. The target application is lymphoid tissue which requires to describe 
lymphocytes. However, the mode of migration is different for these cells. It depends on the so-called 
constriction ring [37-41]. Therefore the active force generation model of fast migrating cells is based 
on the constriction ring model thus taking into account the lateral forces during active migration of 
lymphocytes. 

Cells use a ring to attach a part of their membrane and cytoskeleton to the extra-cellular matrix 
and exert outward directed pressure to their environment by contracting the rear of the cell (Fig. [3]) . 
The ring remains fixed with respect to the extracellular matrix and therefor moves towards the end 
of a cell during migration. In a simplified approach a new ring is generated at the front of a cell only 
when the preceding ring has reached the end of the cell. The force acting on cell i by exerting active 
forces on a neighbor cell j reads 

Fff = a ijP * sign[(x*, - x*) • o*] ^"l,, (10) 

ll x ij x i II 

with the cell orientation Oj , cell surface contact point x*- , constriction ring center x* (which generally 
differs from the cell's center), interaction area and the pressure p* actively exerted by cell i. Note 
that p* is a parameter of the model and shall not be confused with the pressure in ([?]) . 
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Figure 3: Active forces of an actively migrating lymphocyte. A realistic force distribution of a migra- 
tion lymphocyte is shown in the left panel with the constriction ring (thick line) separating 
the area from inwards directed forces from the area with outside directed forces. A mapping 
on a spherical cell representation is shown on the right. The force center is located in the 
center x* of the constriction ring. The plane of the ring is perpendicular to the symmetry 
axis of the force generation determined by the cell's orientation Oj. In the full tissue model 
the spherical surface is approximated by the polyhedral shape of a Voronoi cell determined 
by the neighbor cells j positions such that the active forces acts at the contact points x*-. 

Additionally a constant active force — F^ * (<fo) is directly exerted on the extra-cellular matrix op- 
posite to the direction of the cell's orientation Oj, i.e. the cell is pushed forward against the matrix by 
the force F^ (</»j) in the direction Oj. The orientation of the cell is parallel to the local chemotactic 
gradient (Sec.|23J. 

The motion of cells can also be characterized using a motility coefficient that is similar to diffusion 
coefficient in Brownian motion. However, it is not just a result of physical contacts with the envi- 
ronment but also involves intrinsic features of cell migration. The motility coefficient D is influenced 
by the cell's speed v and the mean free path length [42]. In cells the latter one is determined by the 
persistence time T p during which the cell maintains its polarized structured fixing the direction of 
force exertion and therefore the direction of motion [26-30] . 

D = l -v 2 T p (11) 

Thus the persistence time is a parameter describing the migration of cells which is required in addition 
to the force generation. 

2.3.4. Friction forces 

The component of the friction forces between two cells i and j in contact read 

F- s = jij [xjj — &ij (eij • Xjj)] (12) 

taking into account only the tangential parts of the relative velocity consistent with the treatment of 
the cell-cell interaction by the energy conserving JKR-model ([5]). The force depends on the relative cell 
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velocity Xjj = Xj — Xj. The overdamped approximation and low Reynolds- numbers of cell migration 
justify the linear velocity-dependence [2,4,8]. The unit vector is the normal vector of the plane 
of the cell contact. The friction coefficient 7™ has the dimension of a viscosity times a length scale. 
7ij is chosen proportional to the contact area between the cells. The total drag force is the sum of all 
friction forces F^ rag with neighbor cells plus a term of the free surface of the cell having interaction 
with the surrounding medium. The overall drag force is then given by 



-VmcdRi (l ~ -^j i< (13) 



with medium viscosity r/ me d and the cell-specific viscosities 7ft . Ai = J2jeAf c a *i * s ^ ne sur f a ce m contact 
with other cells. A^ ot = X^'eA/i a ij is the total surface of a cell. The friction coefficient fulfills the 
symmetry 7^ = 7^. The form of the friction coefficients is motivated by the Stokes relation for the 
friction of a sphere at velocity v in a medium with viscosity rj 



pbtokes = Q nr]Rv ( 14 ) 

The geometry related factor 6tt is absorbed in the cell-specific viscosity 7ft in (|13p . The coefficients 
r/ me d and 7ft are chosen such that reasonable active forces have to be generated by the cells in order to 
match the cell velocities measured experimentally (see Table [1]). 



2.4. Reaction-diffusion system of chemoattractants for long-range interactions 

The chemotaxis of cells is described by coupling the direction of migration Oj of a cell i to the local 
chemoattractant gradient. According to the observation that leukocytes tend to have a persistence 
time T p between subsequent orientation changes [20,26-29] the gradient is sensed by the simulated cells 
periodically and Oj is kept constant in between. The concentration of a chemoattractant is computed 
solving the time dependent diffusion equation. The time dependence is included for two reasons. 
First, to account for a dynamics of sources that frequently generates new or removes present sources. 
The time to equilibrate the concentration after changing the source is far longer than the typical 
times for cell migration such that the cells can sense this temporal concentration change. Second, 
the internalization dynamics presented in the next section acts on a time scale comparable to the cell 
migration but much faster than the time scale required for the chemokine diffusion to be equilibrated 
on the tissue scale. 



2.4.1. Internalization model 

The internalization of chemoattractant receptors occurs naturally when they are bound by their 
ligand [43-45]. Studies showed that desensitization of chemoattractant receptors exposed to high 
chemoattractant concentrations [46] occur probably by internalization. Alternative possibilities would 
be a chemical modification of the receptors which is not considered here. Of note, some experiments 
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Figure 4: Receptor desensitization model. The chemoattractant c binds with the rate constant k on 
to its receptor R. The receptor-chemoattractant complex can either dissociate with the rate 
constant k Q R or internalized with the rate constant k\. The internalized complex gets recycled 
re-expressing the internalized receptor R* on the surface. The recycling is characterized by 
the rate constant k r 

fail to detect chemoattractant responses of freshly isolated cells despite the presence of the correspond- 
ing receptor [47]. Thus, the presence of receptors is necessary but not sufficient to cause chemotaxis of 
cells. Most likely the suppression of the receptor function is mediated by a desensitization mechanism 
either by internalization or cross-desensitization [44]. 

Desensitization of cells with respect to a chemoattractant c is achieved by internalization of the 
receptors which have bound the chemoattractant with rate ki [45]. Thus the receptor comes in three 
states: free on the cell membrane (R), bound to the ligand on the cell membrane (Rb), and internalized 
with the ligand (R*) (Fig. 2]), leading to the dynamics: 

R = —k on Rc + k Q f[Rb + k r R* 
Rb = k on R c — k g Rb — hRb 

(15) 

R* = kiR h - k r R* 
6 = —k on R c + k gRb — kc + Q + Ac 

The binding of the ligand is characterized by the on and off rate constants k on and k Q R leading to the 
ligand-receptor complex R^. The rate k v describes the recycling of the internalized receptor R* into 
the free membrane form R. The basic assumption is that the total receptor content i?tot = R+Rb + R* 
is conserved, i.e. no terms describing the transcription or degradation of the receptor are included in 
(fl~5|) . This view is supported by the fact that internalized ligand-receptor complexes dissociate such 
that free receptor become available in the cell [45]. The assumption allows to eliminate the equation 
for R^ from the system (|15p . 

The reaction equations (fT5"j) are completed by an unspecific decay kc of the chemoattractant, a source 
term Q, and a diffusion term Ac for the ligand. The decay term is in accordance with the fact that 
chemoattractants are inactivated and processed by all cells without reaction with the proper receptor 
[44]. This limits the life time of a chemoattractant molecule. The receptors are also transported by 
cells such that the whole system (|15p is a reaction-diffusion system. It is solved using a splitting 
method of first order in time [48,49]. 
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Figure 5: Interaction network regulating migration in the PLF system. The chemoattractants 
CXCL13 and CCL21 guide the migration of B and T cells via the receptors CXCR5 and 
CCR7, respectively. The chemoattractants are produced by different non-migrating cells. 
Thereby FRC can become FDC when the are excited by a B-cell derived signal (LTai/^)- 
FDC return to the FRC state when the signal is absent. 

3. Model of primary lymphoid follicle formation 

This section briefly describes the underlying biological concepts of PLF dynamics. The required 
ingredients to study the formation of the PLF are the B and T cell flow and the dynamics of the 
sessile cell populations, i.e the behavior of FDC and FRC. The parameters for the simulation are 
given in table [TJ A detailed parameter estimate for the reaction-diffusion system (fT5j) is given in the 
appendix. 

3.1. Flux of B and T cells 

B cells are constantly entering SLT predominantly via specialized blood vessels (reviewed in [21]). It 
is less clear how B and T cells leave this tissue. Recent experimental data suggests that lymphatic 
vessels guide cells to leave SLT [21]. In a simplified model cells enter and leave the tissue through 
separate spherical areas that do not further interact with the system. However, one should be aware 
that in real SLT the lymphatic endothelium which may be involved in the efflux of cells is located 
exclusively outside the PLF [50-60]. Thus, there may be some interaction that couples PLF formation 
with the location of exit spots which we ignore in the presented toy model. 

In order to allow the cells to remain and interact within the tissue for a certain time they are not 
permitted to use the exit areas before a minimal transit time has past. This time is about 3 hours 
according to a receptor dynamics [61]. This is consistent with measurements of the minimal transit 
time of lymphocytes [52,62,63]. In order to establish a sufficient cell efflux a chemoattraction towards 
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the exit spots is assumed. This is supported by experimental evidence [21,61,64]. Due to a lack of 
information the corresponding receptor is treated with a simplified model without taking into account 
receptor internalization. 

3.2. Migration of B and T cells 

The migration of B and T cells within SLT follows similar principles. Both cells can perform chemotaxis 
and random motion. It has been shown experimentally that B cells are predominantly attracted by 
the chemoattractant CXCL13 [21]. Both, B cells and T cells respond to CCL21 with T cells being 
more sensitive [21]. The corresponding receptors are CXCR5 (for CXCL13) and CCR7 (for CCL21). 
FRC are the source for CCL21 and FDC produce CXCL13 (FIG.©. 

3.3. Dynamics of sessile cell populations 

Although the origin of FDC is not clear there is strong evidence that they are derived from FRC 
found in the T zone of SLT [65]. Therefore, it is assumed that FDC are an excited state of FRC. The 
excitation is mediated by the contact signal LTai/32 provided by the B cells [12-14] (FIG. [5]). The 
FDC-B cell interaction has been shown to be reversible [12] in the sense that the absence of the B 
cell signal lets the FDC vanish. In the model context that is described as the decay of the excited 
FDC state back into the FRC state. Note that, alternatively to the assumed scenario FRC and FDC 
may also have a common progenitor which can develop in either FRC or FDC depending on external 
stimuli. This is not considered explicitely in the present model. 

Within the PLF model an FRC differentiates to an FDC when the signal threshold for LT«i/?2 has 
been exceeded for a given time Tfrc^fdc- The signal is determined by summing up all LTai/32 con- 
tributions from neighbor cells, i.e. surface density of LTai/?2 times contact area. The differentiation is 
then instantly performed changing the internal cell states of an FRC into that of an FDC, i.e. replacing 
a source for CCL21 by a source for CXCL13. In a similar manner FDC differentiate back to FRC after 
the LTai/32 signal is below the threshold for a critical time Tfdc<— fro F° r simplicity and considering 
the lack of corresponding experimental data Tfrc^fdc = ^FDC<-FRC and equal LTai/32 thresholds 
are assumed. 

3.4. Sequence of events for follicle formation 

The overall picture of PLF formation is as follows. A background of immobile FRC produces CCL21. 
B and T cells enter the tissue through the blood vessels which are represented by few (<10) spheres 
of 30 fan diameter [66] randomly scattered in an area of 150 /xm size. When sufficiently large B cell 
aggregates form, FRC are induced to become FDC by LTai/?2 thus replacing a production site for 
CCL21 by a CXCL13 source. This attracts more B cells and enlarges the forming PLF. T cells are kept 
outside the PLF just by coupling to the CCL21 which is produced around the PLF by the remaining 
FRC. Both, B and T cells, leave the tissue through exit spots which are represented as spheres acting 
as sinks for cells. Any dynamics on the exit spots that may result from PLF formation are ignored 
(compare Sec. 13. 1 [) . 
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parameter 


value 


remarks / ref . 


B/T cell diameter 


9 fim 


[51,67-69] 


Bk 


1 kPa 


[34, 70-72] 




0.4 


[73,74] 


(Tij 


0-0.3 nN^ur 1 


[75,76] 


T P 


120-180 s 


[26,28] 




500 nN /mi" 2 s 


[26-29,70,71,77] 


^FRC— >FDC 


3 h 


[51,78,79] 


D 


10-100 /im' s- 1 


[80] 


size of diffusion grid 


1200 fim 




grid resolution 


35 /im 




max. cell displacement Aa: 


0.9 /im 




min. time resolution At 


10 s 




B : T cell ratio 


0.4 : 0.6 


[63,81] 


influx of cells (B + T) 


0.1-2 cells s -1 


[63,81] 


size of simulation area 


600 /im 




number of FRC 


2500 





Table 1: The parameters used in the simulation for the PLF. The references given support the used 
value. If no references or comments are given, the parameter is a systemic model parameter. 
The value is chosen to guarantee sufficient accuracy of the simulation. The parameters for 
the reaction part of the system (I15D are estimated in detail in the appendix. 

4. Results 

4.1. Biomechanics of lymphoid tissue 

The active forces exerted by lymphocytes to the ECM (F? ct (0,)) and to other cells (p*, (fl"0|) ) are 
not known explicitely (there are just estimates based on experiments of other cell types [35,82-84]). 
The same applies to the friction parameters {rji and ?7med? 

(fT3l ). For simplicity it is assumed that 
Vi = Vmed f° r an lymphocytes. As the known average speed is essentially determined by the ratio of 
active forces and the friction parameters only 2 parameters per cell type remain unknown. These can 
be determined from the measured speed distributions for B and T cells within SLT [26,29]. A reduced 
cell system with a constant number of lymphocytes that respond to a chemotactic point source is 
used. The cell constituents are either B or T cells only or an even mix of both cell types. 

The parameters for the JKR-forces are relatively well known. The parameters entering ([5]) are 
the elasticity of lymphocytes Ei = lkPa [34, 70-72], the Poisson number V{ = 0.4 [73, 74], and 

the surface energy =0 O.SnN/xm" 1 [75,76]. Thus the elastic cell interactions can serve 

as reference forces limiting the reasonable range of friction parameters. Therefore the fit to the 
speed distribution was performed with two qualitatively different friction regimes. In the low friction 
regime [r\i = r] mc d = 50 nN /im -2 s) the active forces that are required to produce the correct speed 
distributions are of the same order as the JKR forces. In this regime it turns out that the mode of 
measuring the cell speed becomes crucial: The speed distribution of the real system was measured 
on the basis of the average displacement over an interval of AT = 10 — 15 s [26,29]. These data 
are reproduced and compared to the simulation results in the two upper panels of FIG. El Using a 
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Figure 6: Comparison of experimental and simulated cell speed distributions. Each bin represents 
an 2 /immin" 1 speed interval. For clear visibility the bins are placed around the centers of 
the speed intervals. The top left panel shows the B cell speed distribution at low (-F^eii = 
18 ± 3nN, p* cell = 0.04 ± O.OlnN^m^ 2 ) and high friction (Fjg^ = 120 ± 20 nN, p*^ = 
1.7 ± 0.05nN/iin~ 2 ) compared to two experimental data sets. The top right panel shows 
the speed distribution of T cells. The way of determining the speed distribution is crucial in 
the case of the low friction (-P^eii = 22 ± 3nN > 2>TceU = °- 06 ± 0.02 nN^m" 2 ). With a high 
sampling rate (AT = 1 s) the distribution has a contribution at high speeds (> 30 /immin -1 ) 
that is not seen using a larger sampling interval (AT = 15 s). For comparison the high friction 
regime (Trf£* u = 240 ± 35 nN, p^ceH = ^-7 ± l.lnN/xm -2 ) is also shown (in that case the 
speed distributions with AT = Is and AT = 15s are virtually indistinguishable). When 
the inter-cellular active forces are omitted (p^cel! = ^Bceii = ([TO]) ) the speed distributions 
show rather sharp peaks (lower panel, sampling interval AT = 15 s). 
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faster sampling with AT = 1 s the speed distribution becomes different and exhibits the influence of 
JKR forces that lead to fluctuations at rather high speed (FIG. [6] top right panel). The JKR forces 
dominantly lead to reversed or perpendicular motion of cells with respect to their migration axis when 
cells collide due to active movement. Thus the net-displacement by JKR forces during AT = 15 s is 
rather small and does only show up at lower speeds in the resulting distribution when longer sampling 
intervals are used. 

In the second regime with high friction (rji = r? me d = 500 nN /im~ 2 s) (FIG. [6]) the average speed with 
AT = Is and AT = 15s match each other (not shown). The resulting speed distribution usually fit 
better to the experimental ones, in particular, because low speed contributions now are shifted to the 
peak resulting in a sharper speed distribution. Moreover the components of the speed distribution in 
the low speed range are smaller as the active forces can more easily overcome the JKR forces such that 
cells are less likely to be immobilized in an cellular environment of high density. However, the active 
forces are higher than measured experimentally [35,82,84]. Yet, the parameter set is computationally 
favorable as the high speed fluctuations are missing. This allows for larger time steps in the simulation. 

The effect of active forces F-j resulting from the constriction ring model is relevant to the speed 
distribution. Without these active inter-cellular forces the speed distribution shows only sharp peaks 
(FIG. lower panel). Also when the speed distribution of each lymphocyte subset is fitted using a 
homogeneous cell aggregate the peaks in the speed distribution tend to be sharper than in experiment. 
This could be compensated with larger active cell interaction forces which however produces the 
wrong speed distribution when B and T cells interact which each other (not shown). 

4.2. Stable follicle formation 

The simple model of PLF formation is able to generate stable follicle sizes of a few 100 /im diameter 
with roughly 10 4 B cells. The follicle forms around the exit spots engulfing them almost completely 
(FIG. [7J (a)). The adjacent T zone is crescent-shaped but tends to form a closed shell around the B 
cells. The T zone is basically determined as non-B zone and T cells occupy the remaining space by 
random migration around the follicle. However, they do not diffuse freely in the whole space due to 
the chemoattraction to the exit spots. Thus this toy model of PLF formation generates a tissue in flow 
equilibrium with dynamically generated sources of cell attraction and keeps B and T cells separated. 

4.3. Effect of receptor internalization 

When studying a stable follicle including chemotaxis towards exit spots the effect of CXCR5 inter- 
nalization dynamics on the follicle shape is not visible. In order to isolate the structural effect of 
CXCR5 receptor internalization, the chemoattraction to the exit spots is switched off. Then a flow 
equilibrium of B cells is no longer possible as the efflux of cells is strongly suppressed (not shown). 
This reflects that randomly migrating cells rarely find the exit spots. Consequently, cell influx and 
efflux are completely shut off which imposes a constant B cell number. This allows to study the 
FRC-FDC dynamics together with the internalization dynamics of CXCR5 alone. It is found that 
the internalization dynamics on its own is a source of instability for the system as the forming follicle 
cannot maintain its shape (FIG. [7] (b)). 

Note, that this result is not affected by considering the internalization dynamics under B cell flow 
conditions. This can be investigated by compensating the low efflux of B cells, which arises due to 
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Figure 7: (Color online) Three-dimensional slice projection of the simulated tissue (a,b). The images 



have been rendered with POVRay ( [http://www.povray.org ). FRC are shown as medium 
sized dark grey spheres (olive), FDC as bright medium size spheres (yellow), and T cells as 
small dark spheres (blue). B cells are shown as small grey spheres (grey to red) (a,b). Cells 
can enter the simulated tissue through blood vessels shown as large grey (red) spheres and 
exit through exit spots indicated as large black (black) spheres. Colors given in brackets 
refer to the colored online version of this figure. 

The localization of the follicle is stable (a). Even a preformed FDC network with a number 
of unsorted B and T cells cannot prevent the formation of the follicle around the exit spots. 
The remainder of the initial configuration is an orientation of the PLF in the direction of 
the initial FDC (a, 8.3 h) which is dissolved later on (a, 12 h). 

The internalization dynamics is visualized in a simulation with a constant number of B 
cells (b). Receptor internalization is found to be a source of instability involving rapid 
morphological changes of B cell distributions. Note that the CXCL13 distribution behaves 
accordingly (c, CXCL13 concentration increases from dark grey (blue) to white (red)). The 
concentration Was visualized with OpenDX ( http ://www.opendx.orgl ). 



switched off chemoattraction to the exit spots, by distribution of the exit spots over a large area. The 
instability is even stronger because the outflowing B cells act as an additional sink for CXCL13 by 
taking surface bound CXCL13 out of the tissue. 

The effect of receptor internalization can be understood by comparing the results to a simulation 
without CXCR5 internalization. If the binding state of the chemoattractant receptors is frozen the 
cells would migrate like being in a potential attracting them to FDC. This is rather independent of the 
B cell-induced dynamics between FDC and FRC because the migration is much faster such that the 
B cell distribution can equilibrate. Thus the distribution of B cells and FDC would match each other 
and the follicle would be stable. The internalization dynamics can however change the 'potential' 
formed by the chemoattractant on a time scale comparable to the cell migration (FIG. 0c)). As the 
cells act as sinks for the chemoattractant they are always generating local gradients away from their 
current position. This cannot be overcome by diffusion of the chemoattractant which is slower than 
the cell migration on the corresponding cellular length scale. Note that the instability only occurs for 
a sufficiently large number of B cells because otherwise the sink for the chemoattractant is too small 
to change CXCL13 gradients significantly. 

The instability of the follicle shape is amplified by the FRC-FDC dynamics. As the B cells spread 
out preferentially at the boundary of the follicle, where CXCL13 concentrations are rather low, they 
extent beyond the follicle border and increase the B cell density outside the follicle. Thus, the area 
covered by the B cells is bigger than expected from a dense packing of the B cells in the FDC area. 
This culminates in the generation of new FDC. As on one hand B cells are very motile and on the other 
hand FDC need some time to dedifferentiate to FRC, on long term this results in an FDC network 
which is bigger than the volume required by the number of B cells. Thus the whole PLF becomes 
unstable in shape following the changed 'center of mass' of the FDC network as determined by the 
concentration peaks of CXCL13. 

5. Discussion 

A model architecture was presented that allows to simulate fast cell migration taking into account 
subtle effects of chemoattractants and cell influx and efflux from tissues. Cells are represented as 
visco-elastic objects. They interact with adjacent cells by passive and active forces. The model allows 
to simulate detailed mechanics of single cells and individual coupling of internal cell dynamics to 
cell mechanics as well as to contact- and long-range-interactions. The use of an underlying regular 
triangulation permits a continuous description of cell positions and motility. It provides an efficient 
way of defining the neighborhood topology for cells of different size in tissue of different density. The 
advantage of the regular triangulation is to provide an effective method to implement fluxes of highly 
motile cells by dynamical modifications of the triangulation. In contrast, lattice-based architectures 
have difficulties to model fast cell migration. Either complicated rules are needed or a lattice-gas model 
is required [85]. In both cases the parameters coupled with the modifications are not easily converted 
to observable quantities. Using a physical representation of cells allows to directly access the involved 
parameters by experiment. This consequently requires a lattice-free description as implemented in 
the present model based on a regular triangulation. The triangulation is complemented by a grid for 
solubles like chemokines. The reaction-diffusion-equations are solved on that grid. Thus, the model 
acquires a hybrid model architecture. The model is based on parameters with physiological well-defined 
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meaning such that the number of parameters is considerably reduced if the corresponding experimental 
values are known. This is, indeed, the case for most parameters in the present test-application. 

The biomechanics of cellular aggregates of fast migrating cells has three components: passive me- 
chanical interaction described by JKR forces, active forces exerted on extra-cellular matrix, and active 
forces directly acting between cells. In general the active forces dominate the JKR forces which allows 
cell movement in dense configurations. The speed distribution with high friction and consequently 
large active forces allows cells to more easily overcome the repulsion generated by JKR forces. Thus, 
the high friction regime mimics large deformations of cells without explicitely taking into account the 
shape of migrating lymphocytes which elongate in the direction of motion [26-30] . Alternatively, the 
elastic modulus Ei of cells could be scaled to appropriate values below the physiologic values and 
correspondingly decreased adhesion parameters <7j. 

The constriction ring model takes into account lateral force generation by actively migrating cells 
which are ignored in other force-based models that consider only the generated net force [2,4]. It is 
shown that the cell interaction due to the exchange of active forces is the major determinant of the 
width of the speed distributions. Thus this interaction might be interpreted as a Brownian motion-like 
fluctuation of the cell speed also it is deterministic. The parameter coupled with the constriction ring 
is the pressure the pressure p* (jlOp . The force i ? i act {(pi) exerted to the extra-cellular matrix may be 
either generated by the integrated pressure of the constriction ring model directly exerted to the matrix 
or by the filopodia dynamics of the classical three-step migration model (reviewed in [77,86]). Recent 
data demonstrates that both migration modes act in lymphocytes [87] although the contribution of the 
constriction ring is greater during chemotactic responses of migrating lymphocytes [87-89]. Especially, 
the characteristic time of the cytoskeleton dynamics of the ring coincides nicely with the persistence 
time [41,87]. The persistence time, the time interval during which the direction of active forces is 
constant, is likely to be an intrinsic feature of the underlying cytoskeleton dynamics [39,41,87,90,91]. 

The test-application for the full model involves a simple approach to the formation of PLF and 
has proven the functionality of the model. More specifically, the motility properties of the cells 
could be adjusted to in vivo motility data of lymphocytes in lymphoid follicles. The model generates 
homeostatic follicles in a flow equilibrium of reasonable size of few 100 fim [50, 92-95] demonstrating 
that the model architecture is able to describe systems that are in a cell flow equilibrium of fast 
migrating cells. 

The effect of CXCR5 receptor internalization was demonstrated under non-flow conditions. Within 
the investigated parameter range (see appendix) the dynamics of CXCR5 and CXCL13 leads to a 
modification of the concentration profile such that cells temporarily migrate away from sources. This 
leads to a frequent change in shape of the cell aggregate. Such a behavior was not observed in vivo [28]. 
Therefore either internalization of CXCR5 does not occur in the PLF system or mechanisms exist that 
counter-balance the shape dynamics. Higher CXCL13 concentrations can overcome this effect as the 
internalization cannot act as a sufficient sink to locally invert the gradient of the chemoattractant 
(not shown). However, the cells then have strongly downregulated receptors such that they loose 
there chemotactic sensitivity, i.e. cells start to migrate randomly preventing efficient aggregation on 
long terms. We observed a stable follicle shape under flow conditions when the B cells chemotax 
to the exit spots without internalization dynamics. Thus, chemotaxis to the exit spots works as 
an independent attractor for B cells that can compenstate for locally inverted gradients of CXCL13 
resulting from the receptor internalization of CXCR5. 

The main purpose of this article was to establish a simulation platform suitable for modeling of 
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homeostatic tissue dynamics on a cellular and subcellular level involving large numbers of fast migrat- 
ing cells. The results have demonstrated that the proposed model design can cope with the complexity 
that occurs in tissue exhibiting a flow equilibrium of fast migrating cells. In particular, it clearly sep- 
arates the various time and length scales and allows to localize the origin of emerging properties on 
the tissue level as well as on the cellular and molecular level. We, therefore, consider this simulation 
tool to be a suitable instrument for the analysis of morphogenesis of highly dynamic tissue. However, 
concerning the application to follicle formation, the assumed mechanisms are only subparts of full 
PLF formation dynamics that need to be refined in order to reproduce microanatomical data and to 
make realistic and quantitative predictions, which is left for future research. 
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A. Estimate parameters for chemoattractant reaction-diffusion system 

The parameters used in the ODE system f|15|) are not known but can be estimated from data of similar 
systems. The dissociation constant for chemoattractants and their receptors are mostly measured 
for other chemoattractants than the ones used here (CCL21 and CXCL13) (Table [2]). The values for 
range from 0.2 nM to 5 nM [96-99]. The dissociation constant for CCR7 with CCL21 has been 
measured to be 1.6 nM [100]. Considering the range of these values it is likely that adopts values 
in this range in the present system. A less favorable situation exists for the reaction rate constant k on 
(k fi = k on K<j). Only few data exists which spread over several orders of magnitude. Values for the 
association rate k on are available for CXCL12 binding to fibronectin [98] (2.5 • 10 5 M _1 s _1 ) . The off 
rate k a s measured for CXCL12 binding to fibronectin is 6.5 • 10 -3 s _1 . 

The two rate constants k\ and k T associated with receptor internalization and recycling (see (|15p ) 
can be estimated from experimental data on receptor desensitization and resensitization experiments 
which are reviewed recently [45]. It is assumed that all non-internalized receptors have bound the 
ligand when high chemoattractant concentrations far above i^d are used as done in the experiment. 
The only equation that is left for the desensitization process is then (|15p (R = 0) 

R* = hiRtot - R*) - k r R* (16) 

with the solution 

R * {t) = TT^ {1 " exp[ " (fci + A;r)t]} - (17) 

The internalized receptor fraction r* at equilibrium becomes 

r * = **(*-°°) = 1 , (18 ) 
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Table 2: Parameter for the equation system (|15|) . 



For the resensitization process we set k\ = and start from = 0) = r*R tot . In the absence of the 
ligand the dynamics for the internalized receptor during resensitization becomes 

R*(t) = r*R tot exp[-k r t}. (19) 

With typical values of r* = 0.3 . . . 0.8 for desensitization and typical recycling times of 60-180 minutes 
(to get r* ~ 0.2 upon resensitization [45]) one arrives at 

k T = 1 • 10" 4 ...7- lO^s- 1 
ki = 5-10 -5 ...3-10~ 2 s _1 . 

Assuming that the internalization process is not in steady state - and solving (|17p numerically - 
doesn't change the results very much compared to the experimental uncertainty. The numerical 
results of PLF- formation are not sensitive to these parameters (data not shown). 

i?tot is not known explicitly. From similar receptors the number of CCR7 molecules on T cells has 
been estimated to be 10 5 per cell [100] and 10 4 for B cells as indicated by the studies that find a 
factor 10 difference of CCR7 levels between B and T cells [101]. Note that the values cited above have 
to be converted into the receptor concentration i? tot based on the cell densities of the corresponding 
lymphocyte type. 
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